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ABSTRACT 

Context. Solar magnetic activity shows both smooth secular changes, such as the Grand Modern Maximum, and quite abrupt drops 
that are denoted as grand minima, such as the Maunder Minimum. Direct numerical simulations (DNS) of convection-driven dynamos 
offer one way of examining the mechanisms behind these events. 

Aims. In this work, we analyze a solution of a solar-like DNS that has been evolved for roughly 80 magnetic cycles of 4.9 years, 
during which epochs of irregular behavior are detected. The emphasis of our analysis is to find physical causes for such behavior. 
Methods. The DNS employed is a semi-global (wedge) magnetoconvection model. For the data analysis tasks we use Ensemble 
Empirical Mode Decomposition and phase dispersion methods, as they are well suited for analyzing cyclic (not periodic) signals. 
Results. A special property of the DNS is the existence of multiple dynamo modes at different depths and latitudes. The dominant 
mode is solar-like (equatorward migration at low and poleward at high latitudes). This mode is accompanied by a higher frequency 
mode near the surface and at low latitudes, showing poleward migration, and a low-frequency mode at the bottom of the convection 
zone. The low-frequency mode is almost purely antisymmetric with respect to the equator, while the dominant mode has strongly 
fluctuating mixed parity. The overall behavior of the dynamo solution is extremely complex exhibiting variable cycle lengths, epochs 
of disturbed and even ceased surface activity, and strong short-term hemispherical asymmetries. Surprisingly, the most prominent 
suppressed surface activity epoch is actually a global magnetic energy maximum, as during it the bottom toroidal magnetic field 
obtains a maximum, demonstrating that the interpretation of grand minima-type events is nontrivial. The hemispherical asymmetries 
are seen only in the magnetic field, while the velocity field exhibits considerably weaker asymmetry. 

Conclusions. We interpret the overall irregular behavior to be due to the interplay of the different dynamo modes showing different 
equatorial symmetries, especially the smoother part of the irregular variations being related to the variations of the mode strengths, 
evolving with different and variable cycle lengths. The abrupt low activity epoch in the dominant dynamo mode near the surface 
is related to a strong maximum of the bottom toroidal field strength, which causes abrupt disturbances especially in the differential 
rotation profile via the suppression of the Reynolds stresses. 
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1. Introduction 

Solar activity manifests itself through the well-known 11-year 
sunspot cycle, during which sunspots are formed progressively 
closer to the equator. This is best seen in the time-latitude do¬ 
main resulting in the so-called butterfly diagram. The cycle is 
not strictly periodic: both cycle length and amplitude are known 
to be variable over time. While there are signs of a long-term 
cyclic modulation on a centennial timescale referred to as the 
Gleissberg cycle, also frequent grand minima with very low 
activity indicators are known in solar history. The Maunder 
minimum (1645-1715) and the Dalton minimum (1790-1830) 
are two prime examples of such minima in existing historical 
sunspot records. Several such events have been retrieved also in¬ 
directly from cosmogenic radionucleid data over several millen¬ 
nia (e.g. |Usoskin et al.|2007| ). Especially the actual duration and 
level of activity of the Maunder minimum (hereafter MM) still 


continues to raise debate. For example, Zolotova_& Ponyavin 


(2015| claimed to have found historical evidence showing that 
some observers did not mark down all the sunspot observations 
on purpose due to the influence of religious or philosophical 
dogmas, resulting in the underestimation of spottedness dur¬ 
ing the MM, which, according to their interpretation, was ac¬ 
tually rather typical cyclic activity during a regular minimum of 
the centennial Gleissberg cycle. In the light of all other avail¬ 
able activity indicators (auroral sightings, cosmogenic radionu¬ 
clide data, solar eclipse observations), analyzed by Usoskin et al. 


( |2015[ ), this interpretation seems however unlikely. Moreover, 
during the MM, the latitude range where sunspots appeared (i.e. 
the width of the butterfly wings) was narrower (e.g. |Ribes & 


Nes me-Ribes| 1 1993 [ |Ivanov &^ Miletsky] |2011| |Usoskin et al. 


201 5|) and strong hemi spherical asymmetry was present (e.g. 


Ribes & Nesme-Ribes[|1993t |Sokoloff & Nesme-Ribes|^994| . 


Analysis of sunspot proper motion seem to indicate slower ro- 
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tation and stronger latitudinal surface differential rotation dur¬ 
ing the MM than for more prominent activity states ([Eddy et al.| 
[19761 |Ribes~& Nesme-Ribes|1993| l. 

The solar magnetic field is thought to arise as an interplay 
of rotation, non-uniformities related to it, and the collective in¬ 
ductive effects of small-scale convective turbulent motions that 
amplify and sustain the magnetic field against intense turbulent 
mixing (see, e.g., Ossendrijver 2003; Charbonneau 2014, and 
references therein). The classical hydromagnetic dynamo picture 
relies on significant turbulence effects throughout the convection 
zone, described by tensorial turbulent transport coefficients de¬ 
scribing, e. g., the a effect, turbulent pumping, and turbulent dif¬ 
fusion (e.g. Moffatt|p~978[ |Krause & Radler|[l980| ). Mean-field 
models of these so-called distributed dynamos usually employ 
subsets of turbulent transport coefficients tha t are either analyt¬ 
ically derived (e.g. [Pipin & Seehafer||2009j) or extracted from 
local convection simulations ( [Kapyla et al.|2006| ). The other dy¬ 
namo paradigm is the so-called flux-transport scenario which 
relies on the existence of highly localized field generation re¬ 
gions. One such region is at the bottom of the convection zone 
in the tachocline, and the other at the top. In the top layer the 
twist of the buoyantly rising flux tubes due to the Coriolis force 
generates poloidal field from the underlying toroidal field (the 
Babcock-Leighton mechanism). The two regions are connected 
through turbulent diffusion and meridional flow acting as a con¬ 
veyor belt (e.g. |Choudhuri et al.|1995[ |Dikpati & Charbonneau 


1999] Karak et al.|2014| ). Both of these scenarios are capable of 

explaining the regular part of the solar cycle, manifested, e.g., 
by the equatorial symmetry properties and the migration direc¬ 
tion of surface magnetic fields, provided the turbulence effects 
are suitably parameterized (e.g. |Charbonneau|2010| ). 

To explain the grand minima-type events with mean-field 
dynamo models, one usually invokes fluctuations in the main 
dynamo drivers, i.e., differential rotation, meridional circula¬ 
tion, and the small-scale turbulence effects. By including such 
fluctuations, dynamo models are usually found to exhibit ir¬ 
regular behavior reminiscent of grand minima-type events (e.g. 


Ossendrij Ver et al.|1996| |Moss et al.|2008| |Choudhuri & Karak| 


2012| . The direct observational knowledge of the change of 
these quantities during a grand minimum-type event, however, 
is very limited, and therefore the mean-field modeling approach 
is problematic. Another possibility is to seek answers from di¬ 
rect numerical simulations of turbulent convection either in local 
(e.g. |Kapyla et al.||2013a|JMasada & Sano||2014|) or global do 


mains (e.g . |Ghizaru et al.|2010( Kapyla et al. 2012; August son 


et al. 2015[|Fan & Fang|2014[|Mabuchi et al.|2015j ; Simitev et al. 


2015| . This is particularly promising in the latter case, where it 
is possible to directly track the change of all relevant dynamo 
drivers, provided that a desired type of dynamo solution is found. 
Important exceptions are the turbulent quantities, namely the in¬ 
ductive and diffusive parts of the mean turbulent electromotive 
force that require special techniques to be properly separated. 
One such method is th e so-called test-field method ( Schr inner| 
et al.|2005l|2007[|MT2| ), the proper application of which to solar¬ 


like glob al, sp herical magnetoconvection solutions is under way 
( [Wamecke et al.|20T6j ). 

Oscillatory dynamo solutions from global magnetoconvec¬ 
tion models have been known for a long time (Gilman 1983; 
Glatzmaier 1985). The first solar-like solutions, however, have 
been obtained only recently ( Kapyla et al. 2012; Augus tson et al.| 
2015), and only a couple of them have previously been run up to 


time scales of interest for detecting the irregular variations: the 
EULAG code Millenn i um simulation at so lar rotation ( |Passos| 
|& Charbonneau||MT4l |Norton et al.||2014[ hereafter EULAG- 


MHD) covers roughly 20 magnetic cycles of 40 years half-cycle 
length, while the ASH code simulation ( Augustson et al. |2015[ 
hereafter ABMT), rotating at three times the solar rate, covers 
roughly 24 cycles with a cycle length of 6.2 years. The for¬ 
mer does not exhibit significant irregular behavior and only very 
weak latitudinal migration of the magnetic field, while the lat¬ 
ter produces a clearer equatorward migrating branch at lower 
latitudes, and a poleward migrating branch at higher latitudes, 
especially pronounced in the radial field. In addition, ABMT re¬ 
port a particular grand minimum-type event, during which mag¬ 
netic activity is suppressed in the equatorial surface region and 
polarity reversals are not seen in the polar surface regions for 
roughly 5 half cycles. The interpretation of the origin of an os¬ 
cillatory magnetic field and equatorward migration vary consid¬ 
erably: from a turbulent dynamo picture producing an aQ os- 
cillatory solution exhibiting a latitudinal dynamo wave ( Kapyla 


et al.||2012} |Warnecke et al. |2014| ) to the magnetic field feed¬ 


ing back to differential rotation via the Lorentz-force, the field 
migration being related to the variations in the differential rota¬ 
tion and possibly to a non-linear dynamo wave ( Augustson et al.| 
[2015] ). 


In this paper we analyze a semi-global (wedge-shaped) 
magnetoconvection simulation similar to those reported earlier 
( [Kapyla et al.|2012[ [Warnecke et al.||2014| ) with slightly varied 
parameters, but integrated over a much longer time span. The 
obtained simulation shows solar-like migration patterns of the 
toroidal magnetic field, and exhibits a dominant cyclic dynamo 
mode with an average cycle length of roughly 5 years. The sim¬ 
ulation covers roughly 80 such cycles. In addition to this ‘ba¬ 
sic’ mode, two other significant modes are detected and char¬ 
acterized with suitable time series analysis tools designed for 
non-periodic signals: Ensemble Empirical Mode Decomposition 
( Wu & Huang|2009| ) and phase dispersion statistics ( |Pelt|1983~| ). 
In the following, we refer to these methods as EEMD and D 2 
statistics, respectively. One of the main goals of this paper is to 
investigate the significance of these multiple dynamo modes to 
the global dynamics of the system. The solution also exhibits 
several epochs of abrupt irregular behavior (disappearance of 
the surface activity, sudden switches of parity, sudden changes 
in cycle length and migration direction of the toroidal field), and 
a physical explanation for these events are sought for by comput¬ 
ing proxies for the different dynamo drivers, namely differential 
rotation, meridional circulation and the a effect, during different 
activity states. 


2. The Model 

Our magnetohydrodynamic (MHD) model has been described in 
many earlier studies, in particular in [Kapyla et al.| ( [2013b| ), and 
the details of it will not be repeated here. We perform compu¬ 
tations in a spherical wedge using spherical polar coordinates 
(r, 0, </>) corresponding to radius, colatitude, and longitude, re¬ 
spectively. The computational domain spans from ro = 0.7 R Q 
to ri = Rq in the radial direction, where Rq = 7 • 10 8 m is 
the solar radius, from Oq = 7r/12 to = ll7r/12 in colati¬ 
tude (±75° latitude), and 90° in longitude, i.e. a quarter of a full 
sphere. Our setup is semi-global in the sense that we exclude the 
polar regions. Including the poles would require prohibitively 
short timesteps. Within this domain, we solve the standard com¬ 
pressible magnetohydrodynamic equations for logarithmic den¬ 
sity In p, specific entropy 5, velocity u , and magnetic vector po¬ 
tential A, which gives the magnetic field as B = V x A. To 
close the system of equations, we assume that the fluid obeys 
the ideal gas law with p = (7 — 1 )pe, where 7 = cp/cy = 5/3 
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is the ratio of specific heats at constant pressure and volume, re¬ 
spectively, and e = cyT is the specific internal energy, where 
T is temperature. The fluid is subject to gravitational accelera¬ 
tion, g = — GM 0 r/r 3 , where G = 6.67 • 10 -11 m 3 kg -1 s -2 
is the gravitational constant and M 0 = 2.0 • 10 30 kg is the 
mass of the Sun, and to rotation, the rotation vector being f2 0 = 
(cos#, — sin#, 0)^o- We neglect self-gravity of the gas in the 
convection zone. 

As explained in detail in Kap yla et al.| ( |2013b| ), since we 
cannot get anywhere near the high Rayleigh numbers of real 
stars, we must use higher diffusivities. In the present model 
this implies a roughly 10 6 times higher luminosity in the model 
in comparison to the Sun. This allows us to reach the Kelvin- 
Helmholtz time scale in our simulations, implying that our runs 
are thermally relaxed. As the convective energy flux scales as 
F CO nv ^ pu 3 , the convective velocity u is roughly 100 times 
greater in the simulations than in the Sun. To obtain the same 
rotational influence on the flow as in the Sun, we must therefore 
increase 9 by the same factor. In general, the scaling of velocity 
and rotation rate can be written as 

Usim = -^ratio^O anC ^ ^sim = -^ratio^O’ W 


where L rat i 0 = Lq/Lq, with Lq and L 0 ~ 3.84 • 10 26 W being 
the luminosities of the model and the Sun, respectively. In the 
current model, as in Warne cke et al.| ( [20T4] ), we then renormalize 
our simulation to the rotation rate 9o = 59 0 , where 9 0 « 
2.7 • 10 -6 s -1 is the mean solar rotation rate, corresponding to 
Q 0 /2 tt = 430 nHz. 

In what follows we express our results in solar units so that, 
say for the velocity, we quote n S i m /L^ t 3 io . The scaling used 
here is based on dimensional arguments. It is supported by mix- 
ing length theory (|Vitense|1953| ) and simulations (Brandenburg 


et al. 2005; Miesc h et al.|2012| ), and should be applicable as long 


as the energy transport is not yet affected by rotation (see e.g. 
Yadav et al.|2013| ). Furthermore, we assume that the density at 


the base of the convection zone (r = O.7f7 0 ) has the solar value 
p = 200 kg m -3 . 

The simulations are performed with the PENCIL COD^j 
which uses a high-order finite difference method for solving the 
compressible equations of magnetohydrodynamics. 


2.1. Initial and boundary conditions 


The initial hydrostatic state is described by a polytrope with in¬ 
dex n a a = 1.5, i.e. an isentropic stratification. Our model does 
not include a stable overshoot region in the bottom of the con¬ 
vection zone. The density stratification is roughly 30 in the be¬ 
ginning, while it decreases to about 20 in the course of the sim¬ 
ulation. In the final state, the number of density scale heights 
covered by the model is roughly 3. To speed up the thermal 
relaxation, the initial condition is chosen not to be in thermo¬ 
static equilibrium, but closer to the final convecting state. We 
choose the heat conduction profile K (r) such that it decreases 
toward the surface like r -15 . This means that a very small por¬ 
tion of the energy is carried by radiative diffusion n ear the sur¬ 
face ( [Brandenburg et al.]|2005t [Kapyla et al.|2014| . We intro¬ 
duce weak small-scale Gaussian distributed noise as velocity and 
magnetic field perturbations in the initial state. 

Our simulations are defined through the energy flux imposed 
at the bottom boundary, F b = —(KdT/dr) | r=ro , the values 
of the rotation rate 9 0 , the kinematic viscosity v , the magnetic 


1 http://github.com/pencil-code 


diffusivity rj , and the subgrid scale diffusivity for the entropy 
in the middle of the convection zone, Xsgs = XSGs(? m = 
0.85 Rq). The radial profile of Xsgs is similar to that shown 
in Fig. 1 of [Kapyla et al.| ( |2011| ) with the surface value being 
Xsgs(^i) = 6xsgS’ corresponding to 6.2 • 10 8 m 2 s -1 in phys¬ 
ical units. The control parameters of the simulation are quan¬ 
tified by the thermal and magnetic Prandtl numbers, PrsGS = 
v/Xsgs an d P r M = v/r), respectively, and the Rayleigh num¬ 
ber Ra = GM 0 (Ar) 4 /cpz/Xg 1 GS R 0 (—d 8 hs/dr) rm , where the 
subscript ‘hs’ indicates a hydrostatic non-convecting state and 
A r = r\ — ro = 0.3 Rq is the depth of the convection zone. 


The radial and latitudinal boundaries are assumed to be im¬ 


penetrable and stress free for the flow; see Eqs. (8)—(9) of Kapyla 
et al. ( |2013b| ). The magnetic field is assumed radial on the outer 


boundary at r i, while the latitudinal and inner radial bound¬ 


aries are perfectly conducting; see Eqs. (10)-(12) of |Kapyld 
[et al.| ( |2013b| ). Density and specific entropy have vanishing first 
derivatives on the latitudinal boundaries. On the outer radial 
boundary we apply a black body condition crT 4 = —K\7 r T — 
Xsgs pT\7 r s, where cr is a modified Stefan-Boltzmann constant, 
which is chosen such that the flux at the surface carries the to¬ 
tal luminosity through the boundary in the initial non-convecting 
state. 


2.2. Diagnostic quantities 


The most important non-dimensional diagnostic quantities 
of our model are the fluid and magnetic Reynolds num¬ 
bers, Re = u rms /uki and ReM = Urms/vh, where 
it rms = yj (3/2 )(v% + v? e ) r Q^t is an estimate of the full three- 
dimensional rms velocity based on the values in the merid¬ 
ional plane, and the subscripts indicate averaging over r, #, 
</>, and a time interval during which the run is thermally re¬ 
laxed, and k{ = 27r/(ri — ro) ~ 21 R^ 1 is an estimate of the 
wavenumber of the largest eddies. The azimuthal velocity com¬ 
ponent is not included in the computation because it is domi¬ 
nated by the differential rotation, which would then not charac¬ 
terize the level of turbulence, which is what we are interested 
in. We also use the meridion al distribution of turbulent veloc¬ 
ities rt' ms (r, #) = yj ( u r2 ) ( f ) t, where the fluctuating velocity is 
defined via u' = u — u, where the overbar denotes an azimuthal 
average, see Sect.[3]for a more specific definition of the averages. 
Furthermore, the rotational influence on the flow is measured by 
the Coriolis number Co = 29o/rt rms &f. 

Our simulations are characterized by Re = ReM ~ 30, 
P r SGS = P r M - 1, Co « 9.5, and Ra « 1.0 • 10 7 . The corre¬ 
sponding numbers computed for the total velocity field, i.e. in¬ 
cluding the azimuthal velocity due to differential rotation, read 
Re to t = R^Mtot ~ 55 and Co to t ~ 5.1. The most impor¬ 
tant input parameters and diagnostic quantities are also listed in 
Tabled] 


This run has therefore a smaller Pr$GS than the runs with 
equ atorward migrating magnetic fields of |Kapyla et al.| ( |2012| ) 
and[Warneck e et al.] ( |2() 1 4[[2() 1 5| ), but twice the value s of Pr$GS 
and PrM than the run with poleward migrating field of Warnecke 


et al.| ( |2014]) by leaving the other parameters the same. The sim- 


ulation of |Augustson et al.| ( |2015| ) has 2-4 times smaller Prandtl 
numbers, but otherwise comparable parameters. The EULAG- 
Millennium simulation setup is similar to that described in 
Ghi zaru et al.| ( |201Q] ) and |Racine et al.| ( |2011| ), but with some ex¬ 
plicit diffusion added for stability reasons. The estimated values 
of their Re and ReM are in the range 30-60, and the magnetic 
Prandtl number is of the order of unity. 
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Table 1 . Input and diagnostic parameters together with the kinetic and magnetic energies realized in the simulation in units of 
10 5 J m -3 . 


PrsGS Pm Re Re M Co Re Mto t E~ E ^ ^ E£g 

1 1 29 29 9.5 55 5.1 4.29 2.45 3.10 1.83 0.96 0.21 0.10 0.66 

Notes. The kinetic energies are defined as the energy of the total flow i?kin = \pu 2 , differential rotation | , meridional circulation 

E^' = 7 ,p(v'r + v'o), and the fluctuating velocity | pu 12 . The magnetic energies are defined as the energy of the total field Smag = 

B 2 /2/xo, azimuthally averaged toroidal E^l g = /3^/2/tn, azimuthally averaged poloidal = (bI+bI) /2/j,o, and the fluctuating component 

C = B j2 /2p 0 . 


3. Data analysis 

Because MHD processes are nonlinear and the resulting data 
non-stationary, it is necessary to choose analysis tools that ac¬ 
curately describe its cyclic components locally and adaptively. 
For example, the solar cycle is known to be of varying length 
and amplitude. While Fourier analysis is the most common data 
analysis technique used to extract periodicities from harmonic 
signals, it requires constant amplitudes and phases and is not 
well-suited to the problem ( jBamhart|2011| ). In this work we have 
chosen to utilize the EEMD and D l statistics, both of which are 
suited for the analysis of non-stationary signals. They are pre¬ 
sented in Appendix |A| and |B| respectively. 

We define mean quantities as azimuthal averages (over the 
^-coordinate) and denote them by overbars. Fluctuations about 
the mean are denoted by a prime, e.g. B' = B — B. We also 
often average the data in time over the period of the simulations 
where thermal energy and differential rotation have reached sta¬ 
tistically saturated states. To study the hemispherical asymme¬ 
tries, we analyze azimuthally averaged quantities separately for 
the northern and southern hemispheres. We also often employ 
the toroidal vs. poloidal decomposition for the azimuthally aver¬ 
aged, and therefore axisymmetric, magnetic field, i.e. write 

B = B tor + i? P oi, (2) 

with B tor = where is the unit vector in the azimuthal 

direction, and where the radial and latitudinal components form 
the poloidal component £? po i = (F> r , £><9, 0). A similar de¬ 
composition is used for the velocity field, u = + u mer , 

where and u mer = (u r ,uo, 0). The three main depths 

at which we plot/analyze the time series are near the surface, 
R s = 0.98 Rq, at the middle, R m = 0.85 i? Q , and close to the 
bottom, i?b = 0.72 R Q , of the convection zone. 


4. Results 


We have run our simulation for nearly 430 years in physical time 
units. The solution we obtain is oscillatory, but it shows very 
complex behavior. The definition of a mean cycle period is not 
unique, but, based on the near-surface toroidal magnetic field os¬ 
cillation, the mean cycle length is, in physical units, 4.9 years, 
and therefore, on average, our simulation covers 80 cycles. In 
comparison to the Sun, the obtained cycle length is roughly by 
a factor of 4 shorter. If, however, this was expressed in solar 
units, the length of our simulation would roughly correspond to 
an evolution of two millennia. From a similar simulation with 
equatorward migrating fields of Warnec ke et al.| p015 ), we find 
that the cycle length increases with slower rotation, showing 
that with more realistic parameters it is possible to obtain cy¬ 
cle lengths comparable to the solar value. To evolve a simu¬ 
lation with a longer cycle over the time span presented here, 


however, is computationally prohibitively expensive. In com¬ 
parison to EULAG-MHD, the number of cycles covered here 
is roughly twice as large, and four times larger in comparison 
to ABMT. Furthermore, more pronounced solar-like latitudinal 
migration, and irregular behavior are observed in our simulation 
than in EULAG-MHD. In the ABMT model, a clear solar-like 
equatorward migration is seen with one MM-type event. Some 
fundamental differences also relate to the data anal ysis, espe- 
cially concerning the EULAG-MHD: In Passos & Charbonneau 
( ]2014[ ) the analogue to the sunspot cycle is defined using the 
toroidal magnetic field strength at the bottom of the convection 
zone, and a poloidal field proxy is built using the radial field 
in the near-surface polar regions, while our analysis covers the 
whole convection zone. Moreover, we restrain ourselves from a 
detailed comparison to the observed statistical properties of the 
solar cycle, as we deem such attempts as somewhat pointless 
due to the vast difference between simulated and real parameter 
regimes, but concentrate on the general effects that may cause 
the irregularities that are seen to occur in our simulation. 


4 . 1. Overall evolution of the magnetic field 

In Fig. [I] we show longitudinal averages of the toroidal magnetic 
field as a function of time and latitude (butterfly diagram) for 
three depths in the convection zone ( R s , R m , and RO, and the 
corresponding plot for the radial field in Fig. [2] The dynamo so¬ 
lution is, in gene ral, very similar to those re ported by Kapyla 
|et al.| ( [2CTT2 ) and |Wamecke et al] ( |2014| |2015| ) with cyclic be¬ 
havior consisting of an equatorward-migrating dynamo wave at 
low latitudes and a polar branch at high latitudes near the sur¬ 
face. The dominating component near the surface is the radial 
one with roughly twice the strength of the azimuthal component, 
while the azimuthal field becomes dominant at larger depths, be¬ 
ing roughly twice (four) times stronger than the radial field in the 
middle ( at the bottom) of the convection zone, see also Figure 
5(e) of Warnecke et al. ( 2014) for a plot of a similar run. The lack 
of strong toroidal magnetic field in the top layer is also related 
to radial field boundary condition at the surface as investigated 
in details in Warneck e et al.| (|2015). The marked difference to 
our earlier results is the irregularity and hemispheric asymme¬ 
try seen in the field evolution. During some epochs, the surface 
activity ceases and/or becomes disturbed (i.e. the cycle length 
is changing significantly), which can happen asynchronously in 
different hemispheres. The radial field evolution is more regular 
than that of the toroidal field, even though the ceased surface ac¬ 
tivity is perhaps even better seen in the radial field evolution, see 
Figs. □ and [2] During the most severe event (20-40 years in the 
south; 35-45 year in the north), both the radial and toroidal fields 
near the surface practically vanish. Simultaneously, regions of 
strong toroidal magnetic field are seen in the near-equator re¬ 
gions (at R s ), and also in the bottom layers (at RO- Even vi- 
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Fig. 1. Time evolution of the mean toroidal magnetic field B$ in the convection zone for three depths (from top to bottom, R s , 
and i?b). The extrema are [—11.0,10.3] kG at R s , [—23.3, 20.7] kG at R m , and [—29.9,31.7] kG at R\>. 


sually, the toroidal field at R^ exhibits another toroidal dynamo 
mode with considerably longer period than the ‘basic’ cycle of 
4.9 years seen throughout the simulation. In addition to these dy¬ 
namo modes, there appears to be a high frequency (short period) 
poleward dynamo mode in the upper part of the convection zone, 
confin ed to the near-equator region, simil ar as in Kapyla et al. 
( |2012| ) and |Warnecke et ak| ( |2014| |2015| ). This dynamo mode 
mig ht be related to a locally operating a z dyna mo in the top lay¬ 
ers ( [Kapyla et al.|2013b]|Warnecke et al.|2014| ). The quantitative 
characterization of these modes is elaborated upon in Sect. |4.3| 


The volume- and time-averaged energies of different flow 
and magnetic field components are listed in Table[l] In Fig.[3ja), 
we show the volume-averaged rms values of the total and mean 
magnetic field together with its mean toroidal and poloidal com¬ 
ponents as functions of time. The dominant component is the 
fluctuating part containing roughly twice as much energy as 
the mean field. In a wedge covering only one quarter of the 
sphere in longitude, the large-scale magnetic field is almost 
purely axisymmetric and contains only a negligible contribution 
from non-axisymmetric modes. Of the mean fields, the domi¬ 
nant component is the toroidal one, containing roughly twice 
the energy of the poloidal component. This is consistent with 
the conclusion of Warnecke et al. ( 2014| that the dynamo is 
of aQ type. The poloidal field evolution shows clear cyclicity 
with modest long-term amplitude variations, whereas the long¬ 


term variations of the toroidal field are much stronger. The ir¬ 
regularities are stronger in the early phases of the simulation, 
the lowered and disturbed surface activity in the butterfly di¬ 
agrams of Figs. [I] and [2] being clearly seen as a global maxi¬ 
mum of the volume-integrated magnetic field energy, especially 
in the toroidal field, see Fig. [3] The early phases, roughly 3/4 of 
the simulation length, show markedly irregular behavior, while 
the latter part, roughly the last 1/4, shows more regular behavior 
with slightly decreased overall magnetic activity level. 


In panels (b) and (c) of Fig. [3j we plot the rms values of 
toroidal (b) and poloidal (c) components for both hemispheres 
(north - blue; south - red). The solid black lines in these pan¬ 
els show the ratio of northern to southern energy (expected to 
attain the value of unity if the hemispheres are totally symmet¬ 
ric). Although averaged over the full time series, the asymmetry 
measure is close to unity for both the toroidal (2.7% dominance 
of the northern hemisphere) and poloidal components (1.7% for 
the north), the variance around the mean is significant (around 
0 .1) for both field components, although the volume-integrated 
quantities indicate no clear correlation between the asymmetry 
and extrema in the total energy. We will investigate the north- 


south-asymmetry in more detail in Sect. 4.6 


In Fig. [3jd) we plot the toroidal magnetic field strength near 
the surface (at R s , red line) and the bottom (at R^, blue line). 
Roughly thrice stronger toroidal fields are located at the bottom 
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Fig. 2. Evolution of the mean radial magnetic field B x in the convection zone for three depths (from top to bottom, R s , R m , and 
i^b). The contours have been clipped at half of the maxima/minima, that is [—11.7,10.4] kG at R s , [—6.1, 5.7] kG at R m , and 
[—3.3,3.8] kG at R^, to make the structures elsewhere than at high latitudes visible. 


of the convection zone in comparison to the surface. The surface 
field shows some variability, but significantly smaller than those 
associated with the bottom toroidal magnetic field. While the 
surface toroidal field strength shows no systematic trend when 
the early and late parts of the simulation are compared, the bot¬ 
tom field is clearly decaying in strength. Therefore, the overall 
decrease and change in the magnitude of the irregularities in the 
total mean field evolution can be attributed to the changes seen 
in the bottom toroidal field. 


4.2. Overall evolution of the velocity field 

In Fig. [4] we show the rms values of the azimuthal, fluctuating 
(a) and meridional (c) velocities. The system energetics is dom¬ 
inated by differential rotation, the energy of which comprises 
roughly 57 percent of the total kinetic energy, being roughly 
2.6 times stronger than the total magnetic energy. The energy 
contained in the meridional flow is negligibly small; in terms 
of the rms velocities, the meridional flow is roughly 15 times 
weaker than the azimuthal one. The energy of the fluctuations 
is roughly 1.3 times weaker than the energy of the axisymmet- 
ric part. Evidently, the azimuthal velocity becomes strongly af¬ 
fected by the dynamically significant magnetic field. A strong 
global magnetic field quenches the azimuthal velocity, while 


during weaker magnetic field epochs, the average azimuthal ve¬ 
locity is higher, as also indicated by the color coding applied 
in Fig. [4^a) and (b) for high, low and nominal states defined 
based on the global magnetic energy level (Dl; see Sect. |4.5| for 
detailed definition). The average meridional flow is an order of 
magnitude weaker, but shows signs of weak systematic variation 
when a running average with the mean magnetic cycle length is 
applied. In Fig.|4][d) we show a zoom-in over a time epoch (50- 
65 years) during which no significant extrema are seen neither 
in the total magnetic field strength nor the azimuthal velocity. 
This figure reveals that both the azimuthal and meridional ve¬ 
locities show a variation over the magnetic cycle, stronger in the 
former, weaker, yet noticeable in the latter. The azimuthal ve¬ 
locity and magnetic field are phase shifted in such a way that 
the azimuthal velocity grows when the magnetic field is weak, 


and decays when it is strong, in broad agreement with Warnecke 
|et al.| ( |2(Tl2] ) and |Karak et al.| ( |2015| ). The meridional velocity 
shows a similar trend, even though it is much more difficult to 
detect as the fluctuations dominate the weak signal. Therefore, in 
terms of the magnetic cycle, the angular velocity variations, of¬ 
ten called torsional oscillations, occur with twice the frequency 
of the magnetic cycle itself (see Fig. [ 7 ] and Sect. |4.3| for a more 
thorough analysis). Our results are in fair agreement with those 
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Fig. 3. (a) Total (T>™ (purple dashed), total mean (T>™, black 

—rms —rms 

solid), toroidal ( B tor , blue), and poloidal (B pol , red) rms values 

as functions of time as averages over the whole volume, (b) B ™ 
in the north (blue) and south (red). The black solid line shows the 
ratio of the northern to southern energy as function of time, and 
the orange line the mean of the ratio computed over the whole 
time span, (c) The same for the poloidal field, (d) The bottom 
(at i?b, blue line) and surface (at R s , red line) toroidal magnetic 
field strengths as functions of time The thin orange horizontal 
lines indicate the mean of the quantities over the whole simu¬ 
lation time. In (a) and (d) thick orange lines show the smoothed 
curve of the total mean magnetic field (a) and the bottom/surface 
toroidal magnetic field strength (d). and the yellow horizontal 
lines the 1.1 and 0.9 times the mean levels, respectively. 




Fig. 4. (a) The rms value of the mean azimuthal velocity as a 
function of time together with the average over the whole time 
series plotted with a black horizontal line. The dashed black line 
shows the volume-averaged rms value of the fluctuating veloc¬ 
ity field, u rms . The red colour indicates the times which we la¬ 
bel as ‘high’, blue symbols ‘low’, and orange ‘nominal’ activity 
states ac cord ing to the strength of the total magnetic field (Dl, 
see Sect. [43}. (b) The rms value of the total magnetic field as a 
function of time, the same colour coding of overplotted symbols 
as in (a). The black horizontal line indicates the mean over the 
whole time series, (c) The rms value of the meridional velocity 
as a function of time (black solid line) with the mean over the 
whole time series (orange horizontal line), and a 4.9-year mov¬ 
ing average (thick red solid line), (d) A zoom-in to 50-65 years 
of evolution of the mean azimuthal velocity (black solid), mean 
meridional velocity (blue line, multiplied by a factor of 10 to 
fit in the figure), and the total mean magnetic field strength (red 
line, scaled to fit the figure). 


therefore not separated for north and south in Fig. [4] It is evident 
that the strong, frequently occurring, short-term asymmetries are 
a property of the magnetic field alone and any asymmetry in the 
flow is probably caused by the asymmetry in the magnetic field. 
In the top panel of Fig. [7] we plot the time-latitude diagram of the 
angular velocity variations at R s . From this plot it can be seen 
that on top of the rather weak, regular oscillation with twice the 
magnetic cycle frequency, which is evident at all depths, there is 
a rather strong long-term variation of the angular velocity near 
the equator in the top layers. The origin of this variation is ana¬ 
lyzed in detail in Sect. |4 5.1| 


obtained for the EULAG-MHD and ABMT models during the 
regular periods of evolution. 

In contrast to the magnetic field, the mean velocity field com¬ 
ponents show only very little hemispherical asymmetry, the evo¬ 
lution of the quantities being virtually identical for both, and are 


4.3. Multiple dynamo modes and their significance 

The results presented in this section have been obtained by ap¬ 
plying EEMD (see App.[A} analysis on the azimuthally averaged 
quantities using an ensemble size of 100 and Gaussian white 
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Table 2. Summary of the mode quantities. 


Period 

Energy fraction 

Mf g 

Tot. Rb Rm Rs 





1 

0.11 

- 

- 

- 

0.13 

0 

0.01 

7 

4.8 

0.46 

0.27 

0.55 

0.39 

0.41 

-0.12 

8 

7.0 

0.11 

0.16 

- 

- 

0.15 

-0.13 

9 

14 

- 

0.11 

- 

- 

0.09 

-0.18 

10 

27 

- 

0.12 

- 

- 

0.18 

-0.38 

11 

53 

- 

0.11 

- 

- 

0.35 

-0.36 

B r 

1 

0.11 

- 

0.17 

- 

- 

0 

0 

7 

5.0 

0.62 

0.29 

0.60 

0.70 

0.75 

-0.10 

11 

52 

- 

- 

- 

- 

0.28 

-0.36 


Be 


1 

0.11 

0.25 

0 

-0.01 

7 

5.0 

0.68 0.44 0.75 0.32 

0.74 

-0.10 

8 

7.4 

0.11 

0.14 

-0.15 

11 

62 

- 

0.37 

-0.58 




6 

2.2 

0.19 

0.25 

0.21 

0.17 

0.03 

0.78 

7 

3.7 

0.14 

0.11 

0.13 

0.15 

0.04 

0.75 

8 

7.2 

0.14 

- 

0.13 

0.14 

0.09 

0.73 

9 

14 

0.11 

- 

- 

0.12 

0.15 

0.70 

10 

26 

- 

- 

- 

0.11 

0.28 

0.74 

11 

48 

- 

- 

- 

- 

0.36 

0.64 


U r 


1 

0.10 

0.61 

0.60 

0.60 

0.62 

0.06 

0.01 

2 

0.19 

0.15 

0.15 

0.15 

0.15 

0.03 

0.02 

11 

41 

- 

- 

- 

- 

0.23 

0.23 


Uq 


1 

0.10 

0.59 

0.60 

0.58 

0.57 

0.03 

-0.05 

2 

0.19 

0.14 

0.14 

0.13 

0.14 

0.01 

0 

11 

43 

- 

- 

- 

- 

0.29 

0.21 


f^kin 


1 

0.10 

0.57 

0.54 

0.56 

0.57 

0.02 

-0.27 

2 

0.19 

0.13 

0.16 

0.13 

0.14 

0.01 

-0.22 

11 

41 

- 

- 

- 

- 

0.38 

-0.48 


Notes. The most prominent modes of mean magnetic and mean veloc¬ 
ity fields found from EEMD. Columns from left to right: mode number, 
average mode period in years, fraction of energy contained in mode 
(over the full meridional plane, over the latitude range at Rb, R m , and 
R s ), regularity level, equatorial symmetry of the modes. To emphasize 
stronger modes we have omitted all energy fraction values lower than a 
preselected level of 10% (indicated by “-”) 


noise with standard deviation equal to that of the time series be¬ 
ing analyzed. The time series were chosen by sampling the simu¬ 
lation domain in radial direction with 12 and in latitudinal direc¬ 
tion with 18 points. Our analysis is focused on the components 
of mean magnetic field, mean velocity field and kinetic helicity, 
aiming to extract the most significant modes for the given time 
series globally as well as at three depths (Rb, R m , and R s ). We 


also use the D 2 statistic (see App.jBj on some selected data sets 
to validate the EEMD results, but also to investigate the coher¬ 
ence of the cycles over time. 

We detected on average a total number of 13 modes. As ex¬ 
pected, it turned out that most of the energy was contained only 
in a limited number of modes. In the following we use a mode 
numbering scheme such that the modes are ordered by their aver¬ 
age frequency from higher to lower values (the highest frequency 
mode thus having an index 1). To express the energy contained 
in each mode quantitatively we will use the notion defined as 
follows: 


r !M)ds 
1 ff(E A f)ds’ 

i 


(3) 


where A is the physical quantity analyzed, i is the mode index, 
A? is the mean square amplitude of the mode, and summation is 
performed over all modes and integration either over latitude at 
selected radii or over the full meridional plane. 

However, calculating only the energy content does not al¬ 
low us to detect weak, but more regular (harmonic) modes. For 
that purpose, for each mode we calculate the spectral entropy 
H =8 Jp(u) log p(y)dv, where p(y) is the power spectral den¬ 
sity at frequency v. Average value of the given quantity per each 
mode compared to the expected value in the case of white noise 
serves as a regularity measure of the mode (i.e. how much more 
Gaussian or non-uniform the spectrum is compared to the spec¬ 
trum of the white noise). An alternative, but more straightfor¬ 
wardly calculated measure is the ratio of energies in each mode 
to the expected ones of white noise. A large value of this ratio 
for a given mode is another indication of the presence of a more 
regular signal. This idea is captured by the following formula: 


m - = E -^r 


(4) 


where Ei is a mode energy fraction defined in Eq. ^ and Ej^ N 
is the expected energy fraction of white noise for the given mode. 
As will be explained later, the results for both regularity mea¬ 
sures are in good agreement with each other. For that reason we 
only show in the following quantitative values calculated based 
on the second definition. 

The results of the EEMD analysis are gathered in Table [2j 
where the mode period represents a descriptive period which is 
calculated based on the zero-crossings of the most dominant in¬ 
stance of the mode of given order (i.e. the intrinsic mode func¬ 
tion, hereafter IMF, with the highest average mean energy over 
the full meridional plane). In the table we also show fractions 
Ei of the energy contained in the it h mode calculated over the 
full meridional plane as well as over the latitude at radii Rb, 
R m and R s . The latter values are approximately equal to the re¬ 
construction error of the initial signal given that the ith mode 
would be omitted from the sum. For clarity reasons only the 
values above 10% are shown, lower values have been marked 
with ith mode would be omitted from the sum. For clar¬ 
ity reasons only the values above 10% are shown, lower values 
have been marked with We can see that in the case of the 
magnetic field components most of the energy is contained in 
mode 7 (hereafter M7), which can be associated with the domi¬ 
nant cyclicity of 4.9 years cycle length deduced from D 2 statis¬ 
tic using surface toroidal field. Only for B$ in the bottom layer 
the energy is more spread between modes with a longer period. 
From Fig. [5j two upper rows showing the spatial distribution of 
the mean amplitudes of the most significant modes found in the 
























































M. J. Kapyla et al.: Multiple dynamo modes and long-term solar activity variations 



Fig. 5. Distribution of mean amplitudes of the most significant modes of meanmagnetic and mean velocity fields found from EEMD. 
The figures are labeled by the physical variable followed by mode index (e.gBg 1 indicates mode 7 of the mean latitudinal magnetic 
field). Colors reflect the mean amplitude (blue - high, yellow - low) of the mode at the given location (we note that the scales of 
separate figures are different). Contours on the plot represent the lines of constant amplitude with values given in the legend. 


magnetic field, it is clear that for the toroidal component of the 
magnetic field, M7 is mainly confined to the region where the 
equatorward migration of the field is observed (compare Fig. [5] 
second panel from left in the upper row with Fig.[l2| top panel). 
For the poloidal field, however, the M7 amplitude is the highest 
at higher latitudes: for the radial field near the surface close to 
the latitudinal boundaries, and for the latitudinal field, near the 
bottom and high latitudes. The modes with shorter cycles (mode 
1 being the one containing most of the energy, hereafter Ml) 


have the highest amplitudes near the equatorial region and near 
the latitudinal boundaries, both regions being concentrated near 
the surface. The long modes (the most notable being mode 11, 
also having an increased regularity measure, hereafter Ml 1) re¬ 
side at the very bottom of the convection zone, confined mainly 
to low-latitude regions. 

In the case of u^, the energy spectrum is more flat with 
roughly equal distribution of energy between modes 6-10. From 
Fig.[5j third row, it is evident that the highest amplitudes of these 
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Fig. 6. Phase dispersion analysis results for the components of B. The calculations were done for at latitude ±22° and radius 
0.94 R q , for B r at latitude ±66° and radius 0.94i? o and for Bq at latitude ±49° and radius 0.82 R Q . Left (right) refers to the 
northern (southern) hemisphere. 


modes all occur in the equatorial region, the region extending 
to higher latitudes for modes 6-7, and getting narrower for the 
ones with a longer cycle. For the other velocity compon ents as 
well as for kinetic helicity, defined as 77kin = c o' u', where 
u/ = V x u', energy is primarily contained only in modes 1 
and 2. From Fig. [5j last row, we can see that these short cy¬ 
cles are most prominent near the surface close to the latitudinal 
boundaries for the radial component of the velocity, and near 
the equator at the surface for the latitudinal velocity component. 
The modes in the kinetic helicity show high amplitudes near the 
equatorial region extending deeper to the convection zone, but 
also in a narrow band near the surface extending the whole lati¬ 
tudinal extent. 

In the penultimate column of Table [2] we show a measure 
of mode regularity as defined in Eq. These values are calcu¬ 
lated at those radii and latitudes where the energy of the given 
mode has its maximum. Higher values of this number indicate 
higher energy content in the given mode compared to the white 
noise level. According to the values shown, besides M7 of the 
magnetic field, we have one additional mode where the regular¬ 
ity has clearly a high value. This is Mil, which, however, has a 
significant energy fraction only in the case of B$ in the bottom 
layer. It is also quite obvious that the leading modes of u r ,uo. 


and i7kin represent primarily noise, because the value of regu¬ 
larity in these cases is very low. It is noteworthy, however, that 
for all these quantities the regularity measure of Ml 1 is high, in¬ 
dicating the presence of a long-term regular cycle. Here we also 
note that the mean spectral entropy calculations lead to similar 
results as above: the spectrum of M7 of the magnetic field is 
the most regular one, while modes 10 and 11 are less regular, al¬ 
though significant deviations from the entropy of the white noise 
spectrum are still detected. It can be concluded, therefore, that 
the dynamo mode at the bottom of the convection zone has an 
influence on the overall dynamics of the system despite it being 
rather insignificant in terms of the total energetics of the system. 

The last column of Table [^represents the equatorial symme¬ 
try of the modes, calculated as a time averaged parity defined by 
Eq. ( p~3] ). The EEMD analysis confirms our earlier conclusion 
that the mean azimuthal velocity field is mostly symmetric be¬ 
tween the hemispheres, as the parity is high. For all the magnetic 
field modes, however, values significantly deviating from perfect 
anti-symmetry (parity values of —1) or perfect symmetry (par¬ 
ity values of +1) are found. The higher frequency modes show 
parity values fluctuating around zero, while the lower frequency 
modes show increasing tendency for antisymmetric parity as the 
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period of the cycle increases. T he e quatorial symmetry will be 
discussed in more detail in Sect. 14.61 

As an independent check for the EEMD results we choose 
the D 2 phase dispersion statistic (see App.[B]). Based on the spa¬ 
tial distribution of the most significant modes, we calculate the 
D 2 spectra at selected depths for the components of B. Ideally, 
we would expect to see a match between the average mode peri¬ 
ods found from EEMD and cycle lengths found from D 2 . 

To focus on the most pronounced activity regions we chose 
the latitude ±22° in the layer at R s for B latitude ±66° in the 
layer at R s for B r , and latitude ±49° in the layer at R m for Bq 
for detecting the short cycles such as M7, see Table [2] Due to the 
limited length of the data set, the D 2 analysis is not feasible for 
detecting the long cycle. However, by using low-pass filtering 
on the raw data, the existence of Mil can be easilyconfirmed. 
The results of the D 2 analysis are depicted in Fig. [6] Evidence 
of M7 with a cycle length of around 5 years can be seen for all 
components. It is interesting, however, that the cycle is modu¬ 
lated differently for northern and southern hemispheres. For the 
southern hemisphere there is an indication of amplitude modu¬ 
lation with a long period, while for the northern hemisphere the 
modulation seems to be more complex. If we assume a simple 
amplitude modulation for the southern hemisphere, and this is 
justified by the fact that we see a relatively symmetric split of 
the spectrum when moving from short coherence times to longer 
ones, then by measuring the difference of these peaks from the 
spectrum taken at the high end of the coherence time, we can 
estimate the period of this modulation. Given the low precision 
that we can achieve using this simple procedure we obtained the 
value for the modulating period to be around 100±10 years. 

We also analyzed the higher frequency region (0.1 to 2 years) 
with the D 2 statistic, but did not detect any strong and/or sin¬ 
gular minima in this range. Taken that the regularity measures 
obtained from EEMD also indicate very low values for these 
modes, we have to conclude that this cyclic mode is coherent 
only over very short time intervals, the variation of the length 
of this cycle being comparable to the length of the cycle itself. 
Therefore, the signal of this cycle can be spread between mul¬ 
tiple EEMD modes, and the spectrum smeared out over a large 
frequency band. 


4 . 4 . Definition of the activity ieveis 


As was evident from the results presented in Sect. |4.1| a de¬ 
creased surface magnetic activity does not straightforwardly im¬ 
ply a global magnetic energy minimum, especially during those 
epochs when the long-period dynamo mode at the bottom of the 
convection zone is strong (roughly speaking the first 3/4 of the 
simulation). Therefore, the definition of a low/high state is not 
obvious. Here we use three different definitions, all based on 
computing running averages with a window coinciding with the 
main periodicity present at different depths and latitudes (4.9 
years): 


Dl. If the volume-averaged mean magnetic energy is 
larger/smaller than the variation of its smoothed aver¬ 
age, the state is termed high/low. All other epochs are 
termed nominal states of activity. 

D2. Otherwise identical, but the surface magnetic field energy at 
R s in comparison to its smoothed value is used as criterion. 

D3. Otherwise identical, but the energy contained in the magnetic 
field mode of the bottom layer at R b is used as criterion. 

In Fig.|3ja; Dl) and (d; D2 and D3) we demonstrate these crite¬ 
ria by overplotting the smoothed curves (orange thick lines) and 


the corresponding upper/lower limits (yellow horizontal lines) 
on top of the original data. In the remaining parts of this section, 
we use these definitions to analyze some key quantities during 
the different activity states. In practise this is done such that we 
compute the mean of a quantity over the whole simulation run, 
divide the data points into different states according to each crite¬ 
rion, and compute the mean profiles for the three activity states 
separately. If the so-obtained profiles differ significantly from 
the global mean, i.e. the difference is larger than the standard 
deviation of the quantity, we take a note of this dependence. 


4 . 5 . Variation of the dynamo numbers 

In this section, we investigate the variability of the three key fac¬ 
tors involved in the dynamo proc ess, na mely differential rotation 
and meridional circulation (Sect. |4.5.1] ), as well as the inductive 
effects due to convective turbulence (a effect; through kinetic 
and magnetic helicities as its proxy, Sect. |4.5.2| ), using the dif¬ 
ferent definitions of the activity level. 


4.5.1. Rotation and its non-uniformities 


As stated in the Sect. [2] the input angular rotation velocity is 
five times the solar rotation rate. The differential rotation profile 
generated in the simulation is solar-like, although the equator is 
rotating only approximately 7% faster than the poles when the 
magnetic fields have grown and saturated at dynamically sig¬ 
nificant strengths. Therefore, the obtained latitudinal differential 
rotation is approximately three times weaker than in the Sun. 
We have also performed a hydrodynamical counterpart simula¬ 
tion, in which roughly twice stronger latitudinal differential ro¬ 
tation (pole-equator difference of 14%) is generated. The dis¬ 
tribution of angular velocity is fairly similar to the MHD state, 
discussed later in this section, and therefore the HD results are 
not shown separately. Furthermore, the hydrodynamic state is 
steady, and shows no oscillatory behavior similar to the dynamo 
run. Therefore, the variations seen in the MHD state (Fig. [4]) 
arise as a consequence of the dynamically significant magnetic 
field to the flow field. Such a backreaction can occur through 
two pathways: a large-scale Lorentz-force J x B (the so-called 
Malkus-Proctor effect; see |Malkus & Proctor119751 ) or through 
turbulence effects either directly on the convective motions and 
thereby the generators of differential rotation (the A-effect: see 
Rudiger 1989) as recently explored in Kara k et ah] (2015|), o r 


through the turbulent Maxwell stress (e.g. |Kapyla et a l. 2004). 
At the same time as quenching the flow field, the magnetic field 
suppresses its own generators. One argument to explain both the 
regular and irregular parts of the solar cycle is related to this 
mechanism (see e.g. ABMT). The classical theory of turbulent 
hydromagnetic dynamos, however, allows for the excitation of 
oscillatory solutions independent of the pre-existence of such 
modulation in the flow field (see e.g. |Parker|1955[ ). In this section 


we set out to investigate the role of the changes seen in the mean 
flow to the long-term evolution of the magnetic field, especially 
to its disturbed states. 


We begin by investigating the angular velocity variations 
over time in more detail. Here S2 V = D — (D) t , where 
Q = Ucft/r sin 6 + £2 0 and (Q) t is time-averaged over the sat¬ 
urated stage. We plot a time-latitude diagram of this quantity at 
R s (Fig. [7] top panel), at which depth the most prominent varia¬ 
tions occur in the equatorial regions. In addition to the strongest, 
seemingly irregular, changes symmetrically distributed around 
the equator, a weaker modulation corresponding to the pole- 
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Fig. 7. From top to bottom: Time-latitude diagram of the angular velocity variations FL, variations of the large-scale Lorentz force 
J x B, and the Reynolds stress components and Rq^ near the surface (at R s ). For the Reynolds stress plots, the contours have 
been clipped at half of the maxima/minima. 


ward migrating magnetic cycle can be seen at high latitudes 
throughout the simulation. In Fig. [7] second panel, we investi¬ 
gate the role of the large-scale Lorentz force J x B in caus¬ 
ing the changes seen in the angular velocity. It is evident that 
the high-latitude variations in the angular velocity, following the 
mean magnetic cycle, are mediated through this effect. The sud¬ 
den drops of the surface magnetic field strength are very clearly 
seen as similar behavior in the large-scale Lorentz force, and 
occur simultaneously in time in comparison to the surface field 
disappearances, especially pronounced during the first 50 years 
of the simulation. Only very mild enhancements of differential 
rotation at higher latitudes are accompanied by the sudden drops 
in the Lorentz force. The strong variations seen in the angular 


velocity near the equator, on the other hand, cannot be explained 
by the large-scale Lorentz force. 

Their symmetric character with respect to the equator hints 
towards them being related to the Reynolds stress component 
R r( f) = u' r u Instead of inspecting all the Reynolds stress com¬ 
ponents, we concentrate here on examining the variability of 
those known to be the most significant for the generation of dif¬ 
ferential rotation, namely R r( p for vertical and Rq^ = u' e u ^ for 
horizontal differential rotation, and postpone the full analysis of 
the turbulent quantities to a forthcoming paper (hereafter Paper 
II). The time-latitude plot of the aforementioned Reynolds stress 
components are shown in the two lowermost panels of Fig. [7] As 
is evident from these figures, both stress components undergo 
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variations on the time scale of the shortest dynamo cycle (Ml), 
while showing no modulation by the dominant magnetic cycle 
(M7). The stresses have non-zero values only near the equatorial 
regions, where the strongest angular velocity variations are seen. 
The minima/maxima of the stresses coincide very well with the 
minima/maxima of the angular velocity and therefore indicate a 
strong connection between these quantities. 

Next we apply the three different definitions of high and low 
magnetic activity states (D1-D3) introduced above, and seek for 
significant deviations, indicative of the sensitivity of the mean 
flow field on the quantity chosen. The mean flow is almost com¬ 
pletely insensitive to the variations of the surface magnetic ac¬ 
tivity level (D2), and not significantly different when the bottom 
activity level is used (D3). Some significant differences, how¬ 
ever, can be detected using the global magnetic activity level 
(Dl) as an indicator. 

Results with Dl are depicted in Fig. [8] where we compare 
the time-averaged rotation profile to the different activity states 
(low, nominal, and high) denoted and defined as 

A^state = ((^)t (^)state)/^0? (5) 

where (VI) t is the time-averaged rotation profile, computed over 
the whole simulation excluding only the initial kinematic stage 
when the magnetic field is still growing, and the quantities in¬ 
cluding angular brackets denote averages over a certain state. 
Using this definition, enhancement/quenching with respect to 
the average value during a certain state will show up as nega¬ 
tive/positive values in the plot, respectively. 

The time-averaged rotation profile is very similar to the ones 
obtained and reported in various earlier works (Kap yla et al.| 
2012, 2013b), being solar-like, but with a somewhat more cylin¬ 
drical profile, showing a mid-latitude region of slower rotation 
leading to a region of reversed sign of radial and latitudinal 
differential rotation. This region was found to be instrumental 
in causing the equatorward migration of the magnetic field in 
| Warnecke et al. | ( |20T4| ) . The variations in the rotation profile av¬ 
eraged over different global magnetic activity states are weak 
(roughly 0.4 percent), while the instantaneous variations near 
the surface (see Fig. [7]) could be as large as 5 percent during the 
early stages of the simulation. During a low/high state, slightly 
faster/slower equatorial rotation is obtained. This is consistent 
with the picture that the stronger the magnetic field, the stronger 
the suppression of the velocity field. The region with a reversed 
gradient of Vt, however, is observed to persist during all activity 
states, with virtually no change in its magnitude. 

In Fig. [8] (middle row, radial differential rotation; lower row, 
latitudinal differential rotation) we investigate how the mag¬ 
nitude and distribution of differential rotation change depend¬ 
ing on the global magnetic activity level. Here we define the 
radial shear as A r = dVt/dr and the latitudinal shear as 
A 0 = dVt/rdO, and the profiles computed over difference ac¬ 
tivity states as 


state — {^i)t (^)z,state? (6) 

where the quantities with angular brackets are averages over the 
whole time series except for the kinematic stage (index t) or 
over a certain state (index i). It is evident that during high/low 
states of magnetic activity, also the radial differential rotation is 
quenched/enhanced similar to the rotation itself, the magnitude 
of the change being roughly 1.5 percent in the equatorial region. 
No significant change is seen for the latitudinal differential ro¬ 
tation, except for at the high latitude regions, which show hemi¬ 
spheric asymmetry especially during the low state. Finally, we 


define a dynamo number describing the magnitude of the radial 
differential rotation as 




A r (Ar) 3 
Vt o 


(7) 


where the turbulent diffusivity is approximated by rjto = 
TU rm S {r,0)/ 3, with r = a M LTH p /u' ims (r,d) being the lo- 
cal turbulent correlation time, c^mlt the mixing length param¬ 
eter and H p = — (<91np(r, 0)/dr)~ 1 the pressure scale height. 
Following our previous studies, we use c^mlt = 5/3. The spa¬ 
tial distribution of Cq follows closely that of the radial gradient 
of the angular velocity (the leftmost panel in the middle row of 
Fig. [8]), and therefore we do not separately plot it, while merely 
report its spatial extrema [—30,110] averaged over the whole 
time span of the simulation, which is similar to those obtained in 
| Kapyla et~ak] ( |2013b| ) and |Warnecke et al.| ( |2014| ). 

Next we examine the meridional flow profiles and compute 
a dynamo number 


a = 


yrrms a 
u mer a/ 


Vto 


( 8 ) 


where u™ s r 


= \/u r 2 + v’l and other quantities are defined in 
the same way as in Eq. 0 - The largest deviation from the time- 
averaged profile is, similar to the differential rotation, obtained 
when the activity states are defined based on the global magnetic 
field energy (Dl), the magnitude of the change being maximally 
6 percent in a thin layer concentrated in the equatorial region 
near the surface. The results are presented in Fig. [9j zoomed 
into the equatorial region. There are, however, some enhanced 
meridional flows generated also near the latitudinal boundaries, 
but as these are most likely artefacts arising from the latitudinal 
boundaries, they are not shown here. The meridional flow pattern 
(Fig.[9]leftmost panel) is very similar to those obtained earlier by 
Kap yld et al.pJl2l|2013b) and |Warnecke' etaL| ( [2QT3][2Ql5| ), i.e. 


consisting of three cells in radius, aligned with the inner tangent 
cylinder. Additionally, there is one anti-clockwise cell confined 
near the surface and equator, the near-surface poleward flow be¬ 
ing the strongest at this location. The strongest variations in mag¬ 
nitude over time are related to this region. The dynamo number 
is varying spatially between 0...3 averaged over the whole sim¬ 
ulation time span, the dynamo number of the meridional flow 
being more than thirty times weaker than of the Vt effect (ra¬ 
dial differential rotation). The high (second panel in Fig. [9]) and 
low (third panel in Fig. [9]) states show no marked difference in 
the spatial distribution nor strength of the main cells, the small 
near-equator cell undergoing the strongest variations. The most 
interesting observation is the marked hemispheric asymmetry 
pronounced especially during the low state: the meridional flow 
in the southern hemisphere is more strongly quenched than the 
northern one. Similar, but the effect being more localized to the 
surface regions, is seen during the high state, when the northern 
surface flow gets reduced, while the southern one gets stronger. 

In conclusion, if the MM was an epoch of a global magnetic 
energy minimum, our results would predict faster surface rota¬ 
tion, stronger differential rotation, and hemispherically asym¬ 
metric meridional flow pattern. If actually a global magnetic en¬ 
ergy maximum was occurring in the deeper layers of the con¬ 
vection zone during the MM, then our prediction would be re¬ 
versed to slower surface rotation and weaker differential rota¬ 
tion, while the meridional flow would retain its asymmetric char¬ 
acter. Neither of these scenarios is consistent with the actual ob¬ 
servations of sunspot proper motion during the MM ( [Eddy et ah 
[1976HRibes & Nesme-Ribes|1993| ). Overall, however, our anal¬ 
ysis suggests that the surface activity disturbances have no direct 
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Fig. 8. Top row, from left to right: normalized time-averaged rotation profile (D) t /D 0 , difference to the rotation profile, AD state , 
during the high, low, and nominal state of global activity (Dl). Middle (bottom) row, the same for mean radial (latitudinal) differ¬ 
ential rotation. 


and/or straightforward relation to changes in the rotational speed 
nor the strength of differential rotation. Moreover, the merid¬ 
ional circulation is so weak that it is most likely incapable of in¬ 
fluencing the global dynamics significantly. During the first 100 
years of the simulation there are strong dips in the angular ve¬ 
locity, accompanied by simultaneous maxima in the global mag¬ 


netic field strength. Only one of these events (during t = 20- 
45 yrs), however, leads to the disappearance of the surface activ¬ 
ity. Even then, the dip in the angular velocity does not coincide 
with the beginning and the duration of the disturbed surface ac¬ 
tivity period. Therefore we can conclude that the surface irreg¬ 
ularities do not originate from the variations seen in the angular 
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Fig. 9. Meridional flow, indicated as arrows, and time-averaged dynamo parameter C u , indicated with color contours. On the left, we 
show the meridional circulation profile averaged over the whole simulation time span, and the next three panels show the difference 
AC Uj (state) = C u — Cu, (state) to the mean during a high, low, and nominal global activity state (Dl). 


velocity or meridional circulation, at the same time noting that a 
disturbed surface activity may not be indicative of a decrease in 
the global magnetic energy level. 


4.5.2. Proxy of the a effect 

We build a proxy for the isotropic a effect following the defini¬ 
tion of |Pouquet et al.| ( [T976] ), 


a = -u' - J' ■ B'/~p), 


(9) 


where the current helicity consists of the fluctuating current 
J r = V x B'/po and the fluctuating magnetic field B' = 
B — B. When looking at the kinetic and magnetic contributions 
to the a effect proxy, it can be observed that the kinetic helicity 
is dominating the signal, the current helicity being roughly an 
order of magnitude weaker. Therefore, we mostly concentrate in 
the following on examining the properties of the kinetic helicity, 
while also computing a suitable dynamo number. 

First we examine the time evolution of the azimuthally aver¬ 
aged kinetic helicity H^ n = c o' u r . This quantity shows a clear 
mean signal as a function of latitude, so that negative (positive) 
values are obtained in the north (south), increasing towards the 
poles where maximal values are obtained, with a strong local¬ 
ized enhancement seen in both hemispheres in the low-latitude 
regions, extending at least half of the depth of the convection 
zone. Close to the bottom of the convection zone at high-latitude 
regions, the sign of this quantity is reversed. These results are in 
agreement with earlier ones both in local and global domains 
(e.g. |Kapyla et al.1|2009l |2013b| ). In Fig. [lOfa), we show this 
quantity with the mean profile over the whole time series sub¬ 
tracted, showing that most of the variation occurs near the equa¬ 
torial region. The signal is dominated by the high-frequency 
modes (Ml and M2), while at higher latitudes the dominant 
mode M7 in the magnetic field is detectable, but not significant 


in comparison to the high-frequency signal. The zoomed-in plots 
of Fig. [TOth) and (c) reveal the presence of the weak basic cy¬ 
cle modulation in the kinetic helicity, but with half of the length 
of the magnetic cycle. Especially during times of clear antisym¬ 
metry, Fig. [l§d), a modulation pattern is seen near the surface 
both at low latitudes (location of the strongest signal in kinetic 
helicity, shown with solid lines in the figure), and at higher lati¬ 
tudes (location of the toroidal magnetic field maximum; dashed 
lines). The modulation pattern is such that the extrema of he¬ 
licity roughly coincide with those of the radial magnetic field 
at mid-latitudes (shown with dotted lines in the figure). During 
times when M7 shows stronger symmetry with respect to the 
equator, Fig. [lOje), all the observed weak dependencies in the 
earlier antisymmetric phase of the simulation break down, al¬ 
though signs of immediate restoration of the modulation pattern 
can be detected when the symmetry abruptly switches at around 
t = 234 yrs back to antisymmetric one. Interestingly, signs of 
a longer-term modulation are seen in Fig. [lOjd), as during that 
epoch the kinetic helicity signal is systematically more negative 
(positive) in the north (south) in comparison to the mean value 
over the whole time span. This is also manifested in Fig. |T0[b) 
as persistent blue (yellow) stripes at low northern (southern) lat¬ 
itudes. Detected by EEMD only through the high value of the 
regularity measure, however, this modulation is of very low en¬ 
ergy in comparison to the high-frequency modes. 

Next, we define and compute a dynamo number 


c a = 


aAr 

Vto 


( 10 ) 


This quantity is plotted in Fig. 11 left panel, as the time average 
over the whole length of the simulation. The profile is similar to 
those obtained earlier, e.g. in Runs Cl and C2 of [Kapyla et al. 
( |2013b| ). 

This dynamo number is almost an order of magnitude larger 
than the corresponding one from meridional circulation, but 
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Fig. 10. (a) Time-latitude plot of the kinetic helicity variations H^ in = H^ n — (H^ n ) t . Zoomed-in plots of H^ in over t = 50- 
65yrs (b) and t = 225-240 yrs (c). (d) at ±10° (solid) and ±25° (dashed) latitude for north (blue) and south (red) with the 
means over the whole time series indicated with an orange horizontal lines. The data is averaged with a running mean over a 1-year 
window, to smooth out the dominant high-frequency component. The dotted lines show the radial magnetic field at ±25° latitude, 
scaled to fit the plot. All quantities are plotted at R s . 


roughly 5 times weaker than the one estimated from radial dif¬ 
ferential rotation. The strongest deviation from the mean profile 
is obtained when high and low states are determined according 
to the surface activity level (D2), while selecting the time points 
based on the global (Dl) and bottom energy (D3) levels produce 
no significant deviation from the mean. During the low state, 
the a effect is strongly enhanced i.e. significantly larger posi¬ 
tive/negative values are obtained in these regions in north and 
south, respectively. During the high state, there is also a mild 
enhancement, while in the nominal state the a effect is strongly 
reduced. The magnitude of the enhancement/quenching is the 
largest of all the dynamo drivers, roughly by 30 percent, but acts 
to the opposite direction as expected, i.e. with a locally and tem¬ 
porally larger C a , the naive expectation is to obtain a more ef¬ 


ficient dynamo near the surface, but weaker surface activity is 
seen. 

We also compute the migration direction of the magnetic 


field predicted by the Parker-Yoshimura sign rule (Yoshimura 
P~976 ) 

^migr = X . (11) 


We plot this quantity on top of the rms-value of the mean toroidal 
magnetic field in the convection zone in Fig.[l2| top panel. In the 
entire region of negative radial shear, manifested by the region 
of slower rotation in the differential rotation profile, equator- 
ward migration is predicted to occur in both hemispheres. This 
region coincides with a toroidal field belt having a maximum 
roughly at ±28° latitude and 0.85 Rq depth. There is another 
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Fig. 11. From left to right: C a time-averaged over saturated stage. Difference to the C a profile during the high, low, and nominal 
state of surface magnetic activity (D2). 


equatorward migration region at the very bottom of the convec¬ 
tion zone, coinciding with the location of the strongest toroidal 
field belt roughly in a similar latitude range. In the upper third 
of the convection zone, however, the predicted migration direc¬ 
tion is poleward. This matches with a toroidal magnetic field belt 
at ±22° latitude and 0.9 i7 0 , depth. No strong radial migration 
is predicted for the lower parts of the convection zone, within 
which both upward and downward regions are seen, but without 
a clear systematic pattern. The good agreement of the predicted 
and actual migration directio n as w ell the mi gration pattern it¬ 
self fit well with the work of |Warnecke et al.| ( 2014[|20l5] ). The 
migration pattern does not change significantly using any of the 
activity level definitions (D1-D3), which is likely due to the fact 
that the gradients of the angular velocity and the kinetic helicity, 
and therefore the a effect, are sensitive to different activity in¬ 
dicators (the former to the global, the latter to the surface) and 
the signal is cancelled out. Therefore, we note that the a effect 
shows the most interesting behavior during the extrema seen in 
the surface activity, but our approach to investigate its effects via 
simple scalar proxies of kinetic and magnetic helicities is inade¬ 
quate; we will return to this problem in Paper II. 

Finally, to investigate whether the different lengths of the 
dynamo solutions and their distinct spatial distribution is due to 
a systematic, strong dependence of the magnetic diffusivity or 
magnetic Reynolds number 


Re M (r, 0) = aMLTg P<msM) , ( 12 ) 

V 

where rj = 10 8 m 2 s -1 is the molecular diffusivity, profile either 
on radius or latitude, we plot their radial profiles at two different 
latitudes in Fig. [12] lower panel. To explain the co-existence of 
dynamo solutions with more than two orders of magnitude vary¬ 
ing periods between the bottom and the top with the dependence 
being due to a weaker diffusivity in the bottom parts, a two or¬ 
ders of magnitude increase of the magnetic Reynolds number as 


function of radius should be seen. Even though the trend seen 
is indeed increasing as function of depth, the magnitude of the 
increase is maximally four, being too weak to provide an expla¬ 
nation. 

4.6. Equatorial symmetry 

In Fig. [13] we plot the equatorial symmetry of the magnetic field, 
defined as 

p _ -Eeven -E'odd ( 13 ) 

^even H - -E'odd 

where F even is the energy of the quadrupolar (symmetric) and 
E’odd the energy of the dipolar (antisymmetric) mode of the mag¬ 
netic field. In Fig.[l3ja) we show the global parity as a function 
of time. This quantity is obtained by computing, for each latitude 
pair for certain depth, the energies contained in the even and odd 
modes, and taking the mean of the obtained data. Evidently, the 
parity is mixed throughout the whole simulation, the dominant 
mode being the dipolar (solar-like) mode (average parity —0.15 
over the whole time series). A difference can again be observed 
between the first 3/4 of the simulation (average parity of —0.26) 
in comparison to the last quarter, when there are larger fluctua¬ 
tions around the clearly positive mean parity of 0.1. 

The best correlation between the different activity level def¬ 
initions can be obtained when D3 is used. The parity not only 
shows a clear modulation with a similar long-term cycle as the 
toroidal field in the bottom of the convection zone, but also a 
clear pattern is revealed for the first 3/4 of the simulation with 
the D3 color coding: during a bottom mode minimum, the par¬ 
ity also obtains a minimum (—0.42), while during a high state, 
the parity is the least solar-like (—0.14). During the last quar¬ 
ter of the simulation, the bottom mode is so weak, that it falls 
below the low activity limit for the whole epoch. Nevertheless, 
the long-term cycle persists both in the bottom mode and in the 
parity, suggesting a relation between these quantities. 
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Fig. 12. Top panel: color contours of B^ in kGauss aver¬ 
aged over the whole simulation run, zoomed-in to show the near 
equatorial region. The overplotted arrows show the predicted 
(Eq. [IT]) time-averaged mean migration direction £ migr normal¬ 
ized to unity. The dashed red lines show the radius at ±25 de¬ 
grees of latitude. Lower panel: radial dependence of the mag¬ 
netic Reynolds number at two latitudes (north - solid; south - 
dashed. 


As can be seen from Table [2] the low-frequency bottom 
mode (Ml 1 ) is predominantly of dipolar symmetry (mean parity 
-0.36). From Fig. [T^b), it is evident that rather strong parity 


fluctuations are related to this mode, the excursions to positive 
parity being more short lived than the periods spent in the nega¬ 
tive parity region. The majority of the parity variation, however, 
is related to the dominant ‘basic’ cycle (M7), which on aver¬ 
age has a parity of roughly — 0 . 12 , but its symmetry is strongly 
fluctuating as a function of time, as is evident from Fig. |T3|b). 
The parity is fluctuating even more strongly than that of Mil, 
between nearly purely dipolar and quadrupolar states. Most of 
the time M7 and Mil and in anti-phase, i.e. when Mil attains 
the largest positive value, M7 is the most negative. There is a ten¬ 
dency for irregular behavior when the parities of these modes are 
in phase, such as during t = 20-45 yrs and t = 250-300 yrs. In 
Fig.[l3jc) and (d) we show two zoom-ins of the parity evolution, 
panel (c) showing a clearly regular, antisymmetric behavior of 

the system and the surface toroidal field (t = 50 -65 yrs), 

while panel (d) shows a predominantly symmetric state that 
during only a half a cycle changes into a antisymmetric one 
(t = 225-235 yrs). During the former epoch, the bottom toroidal 
magnetic field strength attains a minimum, while during the lat¬ 
ter, a maximum is observed. 


4. 7. Variation of general cycle properties 


Evidently, the overall properties of the magnetic cycles show 
quite a large variation over time. In Fig. 14_ we plot three in¬ 
dicators, namely the cycle length of the toroidal field near the 
surface at ±25° latitude, the inclination of the butterfly wing, 
and the butterfly width in degrees for north (triangle symbols) 
and south (pluses) separately. 

As for the toroidal field cycle length (Fig. [14] top panel), two 
distinct types of events can be seen: During suppressed magnetic 
activity near the surface layers (t =20-45, 275-285, and 310- 
320 yrs), the field ceases to change its sign; therefore, some cy¬ 
cles are ‘missed’, and then twice as long cycles are detected as a 
result. The other type of variation are the less abrupt fluctuations 
around the mean cycle length, the shortest cycle lengths being 
slightly less than 4 years and the longest one around around 8 
years. The cycle length is not sensitive to any of the activity level 
definitions (D1-D3), but there is a tendency for the cycles to 
be shorter if the parity is more symmetric (blue symbols), while 
larger cycle length variations occur for more solar-like (antisym¬ 
metric; red symbols) parities. Also, all the ‘missing’ cycles occur 
during the more solar-like parity, during which the bottom cycle 
is attaining its minimum. 

To further analyze the general cycle properties we have fil¬ 
tered the surface magnetic field data so that only points with 
|i? 0 | > 4 kG have been retained. This way the concentrations 
of strong magnetic field are clearly grouped into separate half¬ 
cycles or wings. The relative properties of these structures, such 
as the latitudinal extent (width) and duration can then be directly 
measured. Based on these properties we define the inclination of 
the butterfly wing as the half cycle width divided by its duration. 
For north/south, a solar-like cycle would have negative/positive 
inclinations, while for cycles where migration is weak or non¬ 
existent, values close to zero are expected. As is evident from 
the middle panel of Fig. [14] there is a tendency for the ‘missing’ 
cycles to appear with weak migration only. Anti-solar cycles are 
also observed quite a number of times both for north and south, 
but not at simultaneous epochs. It appears that the shorter the cy¬ 
cle, the higher is the probability for anti-solar migration pattern. 

The corresponding butterfly width (Fig. 14 bottom panel) is 
most often large for the ‘missing’ cycles, which is in striking 
contrast to the observational evidence for the MM (see |Usoskin| 
et al.|2015] and references therein), if these two types of events 
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Fig. 13. (a) Parity P as function of time. The color coding is 
defined using Definition D3. (b) Parity of M7 (black) and Mil 
(red) as function of time; their means are indicated with the hor¬ 
izontal lines with the same color coding, (c) Parity (black) and 
toroidal magnetic field near the surface (blue: north, red: south) 
at ±25° latitude, during the low state of the bottom toroidal 
mode, panel (d) shows the same quantity, but for a high state 
of the bottom toroidal mode. 


are in reality representing one and the same phenomenon. There 
is a tendency for narrow widths to occur when the bottom mag¬ 
netic field is weak. 

4.8. Irregular event during 20-45 years 

Finally, we present a close-up of the most clearly pronounced 
irregular activity epoch during 20-45 years. In Fig. |T5j we re¬ 
produce the plot shown in Fig. [7J but restrict the plotting range 
to 20-65 years, which shows both the irregular evolution, fol¬ 
lowed by the resumption of the regular behavior, indicating that 



J_I_I_I_I_L 
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Fig. 14. Top panel: magnetic cycle length determined from the 
toroidal magnetic field, B at ±25° latitude and r = R s , 
blue/red indicating positive/negative parity, north with triangle 
and south with plus symbols. Middle panel: time evolution of 
butterfly wing inclination. Bottom panel: butterfly wing widths, 
red/orange/blue indicating high/nominal/low bottom magnetic 
field energy (according to D3), north with triangle and south with 
plus symbols. 


the system changes itself into very abruptly after the distortion. 
The time resolution in the figure is high enough to clearly show 
the high-frequency (Ml) mode of the magnetic field (top), not 
well discernible in any other longer time span figure; this cyclic¬ 
ity is seen to persist regardless of the irregularities seen in the 
modes M7 and Mil. 

The irregular behavior is first seen in the southern hemi¬ 
sphere, see Fig.[l5ja), and the event appears to be preceded by a 
field enhancement in the bottom mode (Mil; Fig.[l5{b)). Even 
though the toroidal field of the mode M7 gets very weak, the 
surface activity does not totally cease, but the field evolves with 
significantly longer period without a polarity change. The south¬ 
ern cycle, therefore, appears truly as a ‘missed’ one; meanwhile 
the northern hemisphere continues as usual. After twice the time 
of a normal cycle has elapsed, the bottom magnetic field attains 
another strong maximum. This is reflected by a strong minimum 
in the Reynolds stresses, Fig.[l5jd) and (e), and surface differen¬ 
tial rotation, Fig.p~5]^c). Immediately after this, the toroidal mag¬ 
netic field practically vanishes from the northern hemisphere, 
while the south seems totally unaffected by the northern disrup¬ 
tion. Again, after the south has produced two normal-looking 
magnetic cycles, normal activity is rapidly restored also in the 
north. After that, the bottom magnetic field exhibits no further 
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Fig. 15. From top to bottom: Zoom-ins of the time-latitude diagram of the mean toroidal magnetic field near the surface (at R s ) and 
bottom (at R b), the angular velocity variations, and the Reynolds stress components R r( ^ and Rq$ near the surface (at R s ) during 
t = 20-65 yrs, showing that the magnetic cycle and surface activity are strongly disturbed (t = 20-45 yrs), after which the regular 
state is quickly recovered (t = 45-65 yrs). 


strong extrema, and the mode M7 continues its evolution in a regular manner. It is also noteworthy that the bottom magnetic 

field changes its polarity during the disturbed epoch. 
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During this epoch, it is evident even by eye that there is es¬ 
pecially strong equatorial (north-south) asymmetry related to the 
surface toroidal magnetic field. This is reminiscent of the obser¬ 
vations of the Maunder and Dalton minima ( Ribes & Nesme- 
|Ribes|[T993} |Usoskin et al. |2009] ), that indicate strong hemi¬ 
spheric asymmetry and a relatively strong quadrupolar compo¬ 
nent of the magnetic field. The global parity does not strongly 
reflect this local surface anomaly as a sudden, persistent change 
towards a quadrupolar state (+1), although the quantity fluctu¬ 
ates strongly between the values [-0.9,0.7]. This is most likely 
due to the fact that the symmetry properties of the very strong 
bottom toroidal magnetic field dominate the signal, at all times 
remaining close to a dipolar configuration. 

Although similarly dramatic events cannot be isolated from 
the rest of the time series, it is clear that both extrema and polar¬ 
ity reversals of the bottom mode Mil, frequent during the first 
3/4 of the simulation time, make the system more irregular than 
during the last 1/4 during which especially the energy contained 
in the bottom mode Mil gets weaker. 


5. Conclusions 

In this paper we have analyzed a semi-global (wedge-shaped) 
DNS that produces a solar-like oscillatory dynamo with surface 
migration properties of the magnetic field resembling those from 
observations. The mean cycle length is 4.9 years, and the simu¬ 
lation is evolved over 80 such magnetic cycles. If scaled to solar 
units using the fact that the Sun has a magnetic cycle of roughly 
22 years, our simulation would correspond to two millennia of 
solar evolution. The chosen parameters, most notably the mod¬ 
est stratification and the values of the Reynolds numbers, are 
still far from the real Sun; nevertheless, this combination of pa¬ 
rameters produces a solar-like dynamo, and it is computationally 
affordable to integrate over a large number of cycles. The rota¬ 
tion rate is five times solar, resulting in a roughly three times 
weaker differential rotation than in the Sun. The rotation profile 
is solar-like, but somewhat more cylindrical, and the meridional 
circulation consists of multiple cells in radius. Nevertheless, we 
can expect that this simulation can well be used to study the 
long-term evolution of solar-type dynamos. 

A general property of the dynamo solution is its cyclic na¬ 
ture: even though there is a clear magnetic cycle, the various 
nonlinearities in the system drive it away from an exactly pe¬ 
riodic (harmonic) state. Therefore, special time series analysis 
techniques are needed that can cope with non-periodic signals. In 
this study, we perform a statistical analysis over the whole con¬ 
vection zone, investigating the properties of all relevant quanti¬ 
ties at different depths and latitudes. The methods we choose to 
employ are the EEMD and D 2 statistics. 

The behavior of the dynamo solution is extremely complex. 
In addition to changing cycle length, we observe epochs of dis¬ 
turbed and even ceased surface activity, and strong short-term 
hemispherical asymmetries. The major general findings related 
to the overall dynamics of the system include: 

- The hemispheric asymmetries are related to the magnetic 
field alone, while the velocity field remains almost perfectly 
symmetric at all times. 

- The epochs when the surface activity has decreased or is 
practically non-existent are not global magnetic energy min¬ 
ima but maxima, as strong magnetic fields are stored in the 
deeper parts of the convection zone. 

The main goal of this study is to find the causes for the ir¬ 
regular behavior rather than to make direct comparisons with the 


observed characteristics of the solar magnetic cycle. Therefore, 
we have concentrated our efforts toward analyzing the dynamo 
solution itself to investigate how the most important proper¬ 
ties (cycle length, migration properties, energetics) change dur¬ 
ing the irregular epochs. We also investigated all the key fac¬ 
tors in the dynamo process, namely the rotation and its non¬ 
uniformities, meridional circulation, the inductive action arising 
from turbulent convection (a effect), and the changes of these as 
functions of a relevant activity level measures. The major find¬ 
ings from this part of the analysis can be summarized as follows: 


The specific dynamo solution analyzed here contains not 
only one, but three separate modes, having distinct cycle 
lengths and symmetry properties, and are located in different 
parts of the convection zone. The dominant mode is the near¬ 
surface 4.9 year cycle (denoted with M7) showing equator- 
ward migration at lower latitudes and poleward migration 
at higher latitudes. This mode is accompanied by a weaker 
high-frequency poleward migrating mode (Ml, 0.11 years) 
in the equatorial region and a low-frequency mode at the bot¬ 
tom of the convection zone (Mil, roughly 50 years). 

The crucial property of the different dynamo modes is their 
different symmetries. While the dominating surface mode 
has a mixed equatorial symmetry that undergoes strong fluc¬ 
tuations over time, the bottom mode exhibits nearly pure an¬ 
tisymmetry and is also more regular than the surface mode. 
There is a close relationship between the global magnetic 
and kinetic energies such that strong magnetic fields quench 
the flow field in a manner that the angular velocity is signif¬ 
icantly reduced, differential rotation gets somewhat weaker, 
and the meridional circulation attains an asymmetric char¬ 
acter with modifications of the circulation magnitude espe¬ 
cially in near-equatorial surface regions. We expect that this 
is what would have been seen during the MM, if it was an 
event corresponding to the abrupt disappearance of surface 
activity with a simultaneous global energy maximum in the 
bottom of th e convection zone (see e.g. |Kiiker et al.| fl999; 
Karak 2010, for similar arguments). 


A more common way of interpreting MM is an overall drop 
of the magnetic activity level. Our simulation produces no 
such events, which is naturally no proof of them not occur¬ 
ring in the real Sun. The only quantity markedly linked to the 
surface magnetic activity level in our model is the magnitude 
of the a effect. 

Two kinds of irregularities are separable from the dynamo 
solution itself: smooth variations in the cycle properties and 
abrupt changes that lead to ‘missing’ cycles and ceased sur¬ 
face activity. All these can be satisfactorily explained as the 
interplay of the different dynamo modes and their influence 
on the flow field. 


Even though this simulation has been run for a considerable 
length of time, we still observe a secular trend - especially for 
the long-term cycle in the bottom of the convection zone. To 
fully capture and understand this trend, and to find out whether it 
is a transient or a real cycle, the simulation would still need to be 
continued further. Although we have presented an already rather 
elaborate and time consuming data analysis of the results, we 
have merely touched upon the turbulent quantities. Especially 
the a effect was treated in a very simplified manner, only using 
a proxy from the kinetic helicity, which is far from adequate. We 
recently undertook an effort to compute the full tensorial rep¬ 
resentation of both the diffusive and anti-diffusive contributions 
to the electromotive force using the spherical test-field method 
( Wamecke et al .]|2016| ) from a rather similar, but shorter, solar- 
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like dynamo simulation. It would therefore be useful to compute 
and compare the turbulent transport coefficients from the differ¬ 
ent epochs of the long run presented here. The simplified anal¬ 
ysis presented here, however, re-affirms our earlier conclusion 
of the general dynamo solution and its migration being expli- 
cable in terms of a tu rbulent aVt dynamo (Kapyla et al. 2012| 
|Warnecke et al.|20T4j ). 
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Appendix A: Empirical Mode Decomposition 

Fourier spectral analysis provides a general method for exam¬ 
ining the global energy-frequency distributions. For its validity, 
some crucial conditions must hold: the system must be linear 
and the data must be strictly periodic or stationary. If these re¬ 
strictions are not satisfied, the non-stationary nature of the data 
causes the energy to be spread over a wide frequency range. 
Additionally, deformed wave-profiles, that are a direct conse¬ 
quence of nonlinear effects, need additional harmonic compo¬ 
nents to be fitted. As a result, the energy-frequency distribution 
can be misleading and hard to interpret ( [Huang et al.|1998 ). 

There are different analysis methods developed for the pur¬ 
pose of describing non-stationary signals. Well-known exam- 
ples a re short-time Fourier transform and wavelet transform (e.g. 
Qian 2002; Cohen 1995 ). An alternative to the af orementioned 
methods is the Hilbert-Huang transform (HHT) (Huang et al. 
|1998[ ). In contrast to almost all the other methods, HHT works 
directly in the temporal domain rather than in the corresponding 
frequency space and the basis functions, also known as intrinsic 
mode functions (IMFs), are derived from the data not selected a 
priori. The decomposition makes implicitly the simple assump¬ 
tion that, at any given time, the data may have many coexisting 
simple oscillatory modes of significantly different frequencies, 
one superimposed on the other ( [Huang & Wu|2008] ). 

HHT is performed in two steps: firstly, using an algorithm 
called Empirical Mode Decomposition (EMD), by which the 
signal is separated into a set of IMFs; secondly, for each ex¬ 
tracted mode, Hilbert spectral analysis is applied which allows 
to describe each mode as an analytical signal having the form 
A{t)exp(iip(t )), where A(t) is an instantaneous amplitude and 
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cp(t) an instantaneous phase. For signals that result in such a 
form, the low frequency content is in the amplitude term and the 
high frequency content in the exponential term ( |Cohen|p~995| ). 
Having obtained a decomposition into IMFs which satisfy the 
analytic signal conditions, we can localize any event on the time 
as well as the frequency axis. If local time-dependent aspects of 
the IMFs are not of interest the second step in the above algo¬ 
rithm may be replaced by a crude approach that estimates the 
average mode period and amplitude by counting zero-crossings 
and determining points of extrema. 

One of the major drawbacks of the original EMD was the 
problem known as mode mixing, which is a consequence of sig¬ 
nal intermittency. To overcome this problem, a noise-assisted 
data analysis method called Ensemble EMD (EEMD) was pro¬ 
posed ( |Wu & Huang| |2009). In |Flandrin et al.| ( |2004| ) it was 
shown that when applied to Gaussian white noise, EMD acts 
as a dyadic filter bank. Utilizing this property allows to extract 
robust and statistically significant IMFs. 

EEMD has previously been applied to time series of different 
solar activity proxies. An example of the application to total so¬ 
lar irradiance and sunspot data can be found in |B arnhart] ( |20 11) . 

As an illustration of how EEMD works we have selected a 
time series of B$ taken at R\> and latitude of 22°. This time se¬ 
ries was decomposed into 12 IMFs, most of the energy being dis¬ 
tributed over modes 7-11, each of which individually contribute 
more than 10% to the full energy. The given modes, with the 
corresponding inst antaneous amplitudes are shown on the top 5 
panels in Fig. |A.l| On the bottom panel the original time series 
of Bfj) as well as the sum of these five modes are depicted. The 
small difference between the curves comes from the fact that the 
less significant modes ( 1 - 6 , 12 and 13) have been excluded from 
the sum. We also note that the form of the IMFs is not precisely 
satisfying the criteria of an analytic signal on one hand due to the 
finite precision limit in the IMF extraction algorithm and on the 
other hand due to the relatively small ensemble size (the noise 
does not fully cancel out). 


Appendix B: D 2 phase dispersion statistic 

The D 2 phase dispersion statistic was first introduced in Pelt 
(|1983|). Recent applications of the method can be found in 
ILindborg et al. ( |2013 ), [Karak et al7| ( |2Q15| ), and |01spert et al.| 
( 2015| ). Similarly to EEMD, the D z statistic is well suited for 
non-stationary data. However, while EEMD aims at decompos¬ 
ing the initial time series into a set of quasi-harmonic signals, 
the purpose of D 2 is to find a set of average cycle lengths, P m , 
which are consistent with the full data set. The D 2 statistic is 
defined as follows: 

N—l N 

E 


D 2 (P, At) 


i =1 j=i +1 


E g(ti,tj,P,At)[f(ti) - f(tj)y 


N—l N 

2(j2 E E g(U,tj,P,At) 

i=1 j=i-\-l 


(B.l) 

where f(ti),i = 1,..., N is the input time series, cr 2 is its vari¬ 
ance, g(ti , tj, P, At) is the selection function, which is signifi¬ 
cantly greater than zero only when 

tn ti 


kP , k = ± 1 , ± 2 ,... and 


(B.2) 

|tj -U\ ;$ At, (B.3) 

where P is the trial period and At is the so-called coherence 
time, which is the measure of the width of the sliding time win¬ 
dow wherein the data points are taken into account by the statis¬ 
tic. More precisely, in the given study the function g has been 


chosen as the product of cosine and Gaussian functions. Here 
the former is depending only on the phase difference of given 
points w.r.t given trial period and the latter only on the time dif¬ 
ference between the corresponding points: 

g(ti,tj,P,At) = • g 2 (ti,tj,At), (B.4) 

gi(ti,tj,P ) = 1 ^cos ^2 tt • frac p — ^ + 1^ , (B.5) 

g 2 (ti,tj,At) = exp ln2 (/^~^/~) (B.6) 

As At is made shorter, we match nearby cycles in a pro¬ 
gressively narrower region, and consequently, estimate a certain 
mean cycle length P m , not necessarily coherent for the full time 
span. At large values of At, i.e. when the time window width ap¬ 
proaches the time span of the data, the D 2 spectrum approaches 
the result that would be obtained using e.g. Fourier transform. 
More precisely, P m is determined from the dispersion spectrum 
calculated for the longest coherence time for which the shape of 
the minimum is still symmetric and singular (with a single lo¬ 
cal minimum within reasonably narrow search region to exclude 
possible minima from higher harmonics). Denoting this optimal 
coherence time by At op t we define the mean cycle length as: 

P m = argmin{L> 2 (P, At opt )}. (B.7) 

P 

When visualizing the spectrum, we plot the cycle length P on 
vertical axis, but instead of directly plotting the coherence time 
on horizontal axis, we use At/P instead. This is just a count of 
cycles with a corresponding period that fit into the time inter¬ 
val equal to given coherence time. We call this entity coherence 
length. 
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Fig. A.l. Modes 7-11 of B$ at R b and latitude of 22° 
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